Environmental indicator analysis-Tilefish

# Loading packages and data
library(ggplot2)
library(ecodata)
library(lubridate)
library(dplyr)
library(stringr)
library(marmap) # bathymetry
library(RColorBrewer)
library(ggnewscale)
library(sf)
library(cowplot)
library(tidyverse)
library(ggpubr)
library(sf)
library(ggdist)
library(ggpubr)
library(wesanderson)

#library(glmTMB)

# CPUE data (no env covariates)
gt_data_model_cpue <- read.csv(here::here('data/catch_data/gt_data_model_cpue.csv'))
#names(gt_data_model_cpue) <- tolower(names(gt_data_model_cpue))
# Add in column with cpue 
# note: Paul indicated to use small mesh
gt_data_model_cpue <- gt_data_model_cpue %>% 
       rename_all(., .funs = tolower) %>% 
  mutate(mesh_bin = case_when(mesh_size <= 5.6 ~ 'SM',mesh_size >= 5.6 ~ 'LG',
                                                      TRUE ~ 'NA')) %>%
  mutate(cpue_hr = sum_gt_catch/effort_dur)

# Catch data: 
sfobs <-readRDS(here::here('data/catch_data/gold_tile_sf_ob_v1_temp_price.rds'))

sfob.env <- sfobs %>% 
  mutate(mesh_bin = case_when(mesh_size <= 5.6 ~ 'SM', mesh_size >= 5.6 ~ 'LG',
                              TRUE ~ 'NA'),
         cpue_hr = SUM_GT_CATCH/effort_dur) %>% 
  filter(YEAR %in% c(1998:2022) & mesh_bin == 'SM') %>%
  select(DATE, YEAR, MONTH, YDAY,trip_id,hull_num, area, effort_dur, 
         SUM_GT_CATCH, cpue_hr, mesh_size, depth, start_lat, start_lon, 
         bottomT, bottomT_avg, MIN_TEMP_C, MEAN_TEMP_C, MAX_TEMP_C,
         TEMP_VARIANCE, TEMP_DEVIATION, MEAN_DPTH_M, tri, sed) %>% 
  mutate(YEAR = as.integer(YEAR)) %>% 
  rename_all(., .funs = tolower)


areas <- sort(unique(sfob.env$area))
catch.tally.ann <- sfob.env %>% # aggregate by year
  group_by(year) %>% 
  summarise(ttl_sum = sum(sum_gt_catch))

# Length data from observer program 
lengths <- read.csv(here::here('data/catch_data/gt_data_length_andy.csv')) 
names(lengths) <- tolower(names(lengths))

# Recruitment estimates from 2021 report
recruit <- read.csv(here::here('data/assessment_data/tilefish_rec_estimate_2021.csv'))

# Merge SF/Obs catch data with recruit estimates:
catch_recruit <- cbind(recruit %>% filter(year %in% c(1998:2020)),
                       catch.tally.ann %>%
                         filter(year %in% c(1998:2020)) %>% select(ttl_sum))

# loading in shape files for maps
wd = here::here('shapefiles')
US.areas <- st_read(here::here('shapefiles/USA.shp'), quiet = TRUE)
canada.areas <- st_read(here::here('shapefiles/Canada.shp'), quiet = TRUE)
bts_strata <- st_read(here::here('shapefiles/NES_BOTTOM_TRAWL_STRATA.shp'),
                      quiet = TRUE)
gtf_strata <- bts_strata %>% 
  filter(STRATUMA %in% c('01030', '01040', '01070', '01080', '01110', '01120', 
                         '01140', '01150', '01670', '01680', '01710', '01720', 
                         '01750', '01760'))

bathy <- marmap::getNOAA.bathy(-81,-58, 27, 46)
bathy = fortify.bathy(bathy)

Tilefish data

---

Catch data

Year-class strength is broadly defined as the number of fish spawned or hatched in a given year (Ricker, 1975).

Figure 1. Sum of catch (not accounting for effort), across years. Light blue shaded region represents the temporal range of observer records and red shaded region represents temporal range of study fleet records. The ‘purple’ region is where they overlap. Note that 2000-2005 for observer records had low sample size/number of vessels for tilefish, making the shaded region likely the best region to use for analysis. The vertical dashed lines represent strong year classes for this species (Nesslage et al. 2021). Red asterisk marks year that stock was deemed ‘re-built’.

# tot_catch == total (sum_catch) across hauls. so if tallying up annually, 
# use sum_catch
# Strong year-classes: 1970, 1973, 1993, 1999, 2005, 2013

ggplot(catch.tally.ann, aes(x = factor(year), y = ttl_sum, group = 1))+
   geom_rect(aes(xmin = '2007', xmax = '2022', ymin = -Inf, ymax = Inf), 
            fill = 'red', alpha = 0.02) +
  geom_rect(aes(xmin = '2000', xmax = '2022', ymin = -Inf, ymax = Inf), 
            fill = 'lightblue', alpha = 0.05) +
   geom_vline(xintercept = c('1993','1999', '2005', '2013'), lty = 2) +
  geom_line(color = 'black', size = 1.5) +
  annotate("text", label = "*",
    x = 26, y = 14000, size = 8, colour = "red" )+
  xlab('Year') + 
  ylab('Total sum tilefish catch') + 
  # facet_wrap(~month)+
  theme(axis.text.x = element_text(color = 'black',
                                   size = 12, angle = 45, vjust = 1, hjust=1)) +
  ecodata::theme_facet()

CPUE

Figure 2. Catch-per-unit-effot for undirected trawl trips from the Study fleet and observer program. Zeros have been added using species association methodology (via jaccard index).

see here for example

gt_data_model_cpue %>% 
  filter(mesh_bin == 'SM') %>% # note: Paul indicated to use small mesh
  group_by(year, source) %>% 
  summarise(mean_cpue = mean(cpue_hr),.groups = 'drop') %>%
  ggplot(aes(x=year,y=mean_cpue)) +
  geom_line(lwd = 1) +
  facet_wrap(~source) + 
  theme_bw()

gt_data_model_cpue %>% 
  filter(mesh_bin == 'SM') %>%
  group_by(year) %>% 
  summarise(mean_cpue = mean(cpue_hr),.groups = 'drop') %>%
  ggplot(aes(x=year,y=mean_cpue)) +
  geom_line(lwd = 1) +
  labs(title = 'Study fleet + Observer combined') + 
  theme_bw()

Maps (catch)

Tilefish catch locations (study fleet/observer)

yrs = sort(unique(gt_data_model_cpue$year))    

#for(i in 1:length(yrs)){
yrmap <- function(yrs){
  gt_data_model_cpue %>% 
  filter(start_lat < 42.5 & depth_est > 50 & year == yrs) %>%
  mutate(bin = cut(year, seq(min(year), max(year) + 4, 4), right = FALSE)) %>%
  ggplot() + 
  geom_sf(data = US.areas %>% st_as_sf(),color = 'gray20', fill = '#cbdbcd') +
  geom_contour(data = bathy,
               aes(x=x,y=y,z=-1*z),
               breaks=c(50,100,150,200, Inf),
               size=c(0.3),
               col = 'darkgrey') +
  stat_summary_2d(aes(x=start_lon, y=start_lat, z = cpue_hr),
                  binwidth=c(0.16666,0.16666)) + 
  scale_fill_viridis_c() + 
  theme(legend.position = "bottom",
        legend.key.size = unit(0.2, "cm"),
        legend.key.width = unit(1, "cm")) +
  coord_sf(xlim = c(-75,-65.5), ylim = c(36,44), datum = sf::st_crs(4326))  +
  labs(x = '', y = '', fill = 'CPUE') +
  theme_bw() 
}


for(i in 1:length(yrs)){
 cat("\n####",  as.character(yrs[i]),"\n")
    print(yrmap(yrs[i])) 
    cat("\n")   
}

2000

2001

2002

2003

2004

2005

2006

2007

2008

2009

2010

2011

2012

2013

2014

2015

2016

2017

2018

2019

2020

2021

2022

Recruitment estimate

Figure 3. Age-1 recruitment estimate from the 2021 tilefish assessment across all years

ggplot(recruit, aes(x = factor(year), y = recruit_est, group = 1))+
   geom_rect(aes(xmin = '2007', xmax = '2022', ymin = -Inf, ymax = Inf), 
            fill = 'red', alpha = 0.02) +
  geom_rect(aes(xmin = '2000', xmax = '2022', ymin = -Inf, ymax = Inf), 
            fill = 'lightblue', alpha = 0.05) +
   geom_vline(xintercept = c('1993','1999', '2005', '2013'), lty = 2) +
  geom_line(color = 'black', size = 1.5) +
  annotate("text", label = "*",
    x = 26, y = 14000, size = 8, colour = "red" )+
  xlab('Year') + 
  ylab('Total sum tilefish catch') + 
  # facet_wrap(~month)+
  theme(axis.text.x = element_text(color = 'black',
                                   size = 12, angle = 45, vjust = 1, hjust=1)) +
  ecodata::theme_facet()

### Thought: Should we isolate years associated w/strong year classes (or bad)
###           for correlations and analyses? 

Figure 4. Recruitment estimates in focus years

• Paul has suggested that data used to calculate recruitment estimates was strongest from 2000 to present.
• Extended to 1998 to capture both the strong year class and El nino associated with that year.

ggplot(recruit %>% filter(year %in% c(1998:2022)),
       aes(x = factor(year), y = recruit_est, group = 1))+
  geom_vline(xintercept = c('1993','1999', '2005', '2013'), lty = 2) +
  geom_line(color = 'black', size = 1.2) +
  xlab('Year') + 
  ylab('Recruit estimates') + 
  theme(axis.text.x = element_text(color = 'black',
                                   size = 12, angle = 45, vjust = 1, hjust=1)) +
  ecodata::theme_facet()

Figure 5. Recruitment estimates and Study Fleet and Observer catch data. Black line denotes recruitment estimate, yellow denotes sum of annual catch data across both Study fleet and Observer programs.

options(scipen=999)
ggplot(catch_recruit) +
  geom_line(aes(x = factor(year), y = recruit_est, group = 1),
            col = 'black', size = 1.2) +
  geom_line(aes(x = factor(year), y = ttl_sum*1000), size = 1.2, 
            color = 'goldenrod1', group = 1) +
  scale_y_continuous(sec.axis = sec_axis(~./1000, name = 'Catch (lbs)')) + 
  geom_vline(xintercept = c('1993','1999', '2005', '2013'), lty = 2) +
  xlab('Year') + 
  ylab('Recruit estimates') + 
   theme(axis.text.x = element_text(color = 'black',
                                   size = 12, angle = 45, vjust = 1, hjust=1)) +
  ecodata::theme_facet()

Length data

All lengths

Figure 6. Distribution of lengths Figure 7. Length frequencies Figure 8. Frequency of smaller individuals

# Define category breaks
size_breaks <- c(0,10,20,30,40, 50, 60, 70, 80, 90, 100)
# Making a function to bin the catches
label_interval <- function(breaks) {
  paste0("(", breaks[1:length(breaks) - 1], "-", breaks[2:length(breaks)], ")")
}
labels = label_interval(size_breaks)

# length freq. table
tab = table(cut(lengths$lenanml, 
    breaks = size_breaks,
    labels = label_interval(size_breaks)))

## Plot full distribution
ggplot(lengths,
       aes(x = lenanml)) + 
  geom_bar(position = position_dodge(), 
           alpha = 0.4, fill= 'blue', color="black") + 
  xlab('Tilefish length (cm)') +
  theme_bw() +
  theme_facet()

# Plot length frequencies
barplot(tab, xlab = 'Length bins (mm)', main = '')

# Just the little ones
barplot(tab[1:3], xlab = 'Length bins (mm)', main = '')

Juveniles

Young of year - year 1 and 2 size class

ggplot(lengths %>% filter(lenanml <= 26),
       aes(x = lenanml)) + 
  geom_bar(position = position_dodge(), 
           fill= 'slateblue', color="black") + 
  xlab('Tilefish length (cm)') +
  theme_bw() +
  theme_facet()

ggplot(lengths %>% filter(lenanml <= 26),
       aes(x = lenanml, fill = numlen)) + 
  geom_bar(position = position_dodge(), 
           alpha = 0.4, fill= 'blue', color="black") + 
  xlab('Tilefish length (cm)') +
  theme_bw() +
  facet_wrap(~year) + 
  theme_facet()

Environmental data

---

The strong year classes for Golden Tilefish were 1993, 1998, 2005, 2013. Some of the underlying oceanographic processes that may be related to recruitment may influence habitat, retention/displacement and food availablity. These are explored below.

Habitat

Tilefish occupy a very narrow band of habitat conditions. Therefore, temperature and salinity may be of interest.

SST

# SST

SST Fronts

# SST fronts

Bottom Temperature

Figure 1. GLORYS vs in-situ bottom temperatures from study fleet vessels.

Figure 2. Bottom temperature (C) across years. Blue dots are in-situ data, red dots are from GLORYS.

ggplot2::ggplot(sfob.env, aes(x=bottomt, y=mean_temp_c)) +
  geom_point(color="blue", alpha=0.1)+
  geom_abline(intercept = 0, slope = 1) +
  xlab('Bottom Temp (SF)') +
  ylab('Bottom Temp (GLORYS)') +
  theme_bw() 

ggplot2::ggplot(sfob.env, aes(x=bottomt, y=year)) +
  geom_point(color="blue", alpha=0.1) +
  geom_point(data = sfob.env, aes(x=mean_temp_c, y=year),
             color="red", alpha=0.1) +
  xlab('Bottom Temp') +
  ylab('Year') +
  labs(color = 'Source') +
  theme_bw() 

Maps (Bottom T)

jet.colors <-colorRampPalette(c("blue", "#007FFF", "cyan","#7FFF7F", "yellow", "#FF7F00", "red", "#7F0000"))

# select just years with study fleet bottom temps
sf.bt <- sfob.env %>% filter(year>2006 & depth > 50) 
yrs = sort(unique(sf.bt$year))    

#for(i in 1:length(yrs)){
yrmap <- function(yrs){
  sf.bt %>% filter(year == yrs) %>% 
  ggplot() + 
  geom_sf(data = US.areas %>% st_as_sf(),color = 'gray20', fill = '#cbdbcd') +
  geom_contour(data = bathy,
               aes(x=x,y=y,z=-1*z),
               breaks=c(50,100,150,200, Inf),
               size=c(0.3),
               col = 'darkgrey') +
  stat_summary_2d(aes(x=start_lon, y=start_lat, z = bottomt),
                  binwidth=c(0.16666,0.16666)) + 
    scale_fill_gradientn(colors = jet.colors(20)) +
  coord_sf(xlim = c(-75,-65.5), ylim = c(36,44), datum = sf::st_crs(4326))  +
  labs(x = '', y = '', fill = 'Bottom temperature (°C)') +
  theme_bw() 
}


for(i in 1:length(yrs)){
 cat("\n######",  as.character(yrs[i]),"\n")
    print(yrmap(yrs[i])) 
    cat("\n")   
}

2007

2008

2009

2010

2011

2012

2013

2014

2015

2016

2017

2018

2019

2020

2021

2022

Note here temperatures are averaged across all depths > 50 for each month.

# Note here temperatures are averaged across all depths > 50 for each month.

# Create in-situ bottom temps by month w/lag
df.lag = sfob.env %>% filter(year > 2006 & depth > 50) %>%
  group_by(year,month) %>% 
  summarise(mean_dpth = mean(depth),
            mean_bt = mean(bottomt)) %>% 
  mutate(mean_bt_lag2 = lag(mean_bt,2),
         mean_bt_lag3 = lag(mean_bt,3), 
         mean_bt_lag6 = lag(mean_bt, 6))
# Join in-situ bottom temps w/assessment recruitment estimate
df.join = dplyr::full_join(recruit, df.lag, by = join_by(year)) %>%
  dplyr::select(year, month, recruit_est, mean_dpth, 
                mean_bt, mean_bt_lag2, mean_bt_lag3, mean_bt_lag6) %>%
  tidyr::drop_na()
# See what months have data
sort(unique(df.join$month))

## [1]  7  8  9 10 11 12

hist(df.join$month) # will group into spring/summer fall/winter categories

## spring/summer bottom temp no lag
ggplot2::ggplot(df.join %>% filter(month %in% c(4,5,6,7,8)),
                aes(x=recruit_est, y=mean_bt)) + 
  geom_point(color = 'black') +
  stat_smooth(method = "lm",
              formula = y ~ x,
              geom = "smooth") +
  xlab('Recruitment estimate') +
  ylab('Mean bottom temp (°C)')+
  stat_cor(aes(label=..rr.label..)) +
  theme_bw()

## spring/summer bottom temp no lag
ggplot2::ggplot(df.join %>% filter(month %in% c(9,10,11,12)),
                aes(x=recruit_est, y=mean_bt)) + 
  geom_point(color= 'black')+
  stat_smooth(method = "lm",
              formula = y ~ x,
              geom = "smooth") +
  xlab('Recruitment estimate') +
  ylab('Mean bottom temp (°C)')+
  stat_cor(aes(label=..rr.label..)) +
  theme_bw()

## spring/summer bottom temp 2 month lag
ggplot2::ggplot(df.join %>% filter(month %in% c(4,5,6,7,8)),
                aes(x = recruit_est, y = mean_bt_lag2)) + 
  geom_point(color = 'black') +
  stat_smooth(method = "lm",
              formula = y ~ x,
              geom = "smooth") +
  xlab('Recruitment estimate') +
  ylab('Mean bottom temp (°C)')+
  labs(title = 'Lag 2 months') +
  stat_cor(aes(label=..rr.label..)) +
  theme_bw()

## spring/summer bottom temp 2 month lag
ggplot2::ggplot(df.join %>% filter(month %in% c(9,10,11,12)),
                aes(x = recruit_est, y = mean_bt_lag2)) + 
  geom_point(color= 'black')+
  stat_smooth(method = "lm",
              formula = y ~ x,
              geom = "smooth") +
  xlab('Recruitment estimate') +
  ylab('Mean bottom temp (°C)')+
  labs(title = 'Lag 2 months') +
  stat_cor(aes(label=..rr.label..)) +
  theme_bw()

## spring/summer bottom temp 3 month lag
ggplot2::ggplot(df.join %>% filter(month %in% c(4,5,6,7,8)),
                aes(x = recruit_est, y = mean_bt_lag3)) + 
  geom_point(color = 'black') +
  stat_smooth(method = "lm",
              formula = y ~ x,
              geom = "smooth") +
  xlab('Recruitment estimate') +
  ylab('Mean bottom temp (°C)')+
  labs(title = 'Lag 3 months') +
  stat_cor(aes(label=..rr.label..)) +
  theme_bw()

## spring/summer bottom temp 3 month lag
ggplot2::ggplot(df.join %>% filter(month %in% c(9,10,11,12)),
                aes(x = recruit_est, y = mean_bt_lag3)) + 
  geom_point(color= 'black')+
  stat_smooth(method = "lm",
              formula = y ~ x,
              geom = "smooth") +
  xlab('Recruitment estimate') +
  ylab('Mean bottom temp (°C)')+
  labs(title = 'Lag 3 months') +
  stat_cor(aes(label=..rr.label..)) +
  theme_bw()

## spring/summer bottom temp 6 month lag
ggplot2::ggplot(df.join %>% filter(month %in% c(4,5,6,7,8)),
                aes(x = recruit_est, y = mean_bt_lag6)) + 
  geom_point(color = 'black') +
  stat_smooth(method = "lm",
              formula = y ~ x,
              geom = "smooth") +
  xlab('Recruitment estimate') +
  ylab('Mean bottom temp (°C)')+
  labs(title = 'Lag 6 months') +
  stat_cor(aes(label=..rr.label..)) +
  theme_bw()

## spring/summer bottom temp 6 month lag
ggplot2::ggplot(df.join %>% filter(month %in% c(9,10,11,12)),
                aes(x = recruit_est, y = mean_bt_lag6)) + 
  geom_point(color= 'black')+
  stat_smooth(method = "lm",
              formula = y ~ x,
              geom = "smooth") +
  xlab('Recruitment estimate') +
  ylab('Mean bottom temp (°C)')+
  labs(title = 'Lag 6 months') +
  stat_cor(aes(label=..rr.label..)) +
  theme_bw()

Salinity

Here we explore salinity from the GLORYS reanalysis model at three different depths

• 55 meters (relevant to larvae, recruits)
• 110 meters (relevant recruits, juveniles)
• 220 meters (relevant to juveniles, adults)

# Salinity 

Currents

Cross-shelf processes may influence the retention or displacement of tilefish during early life history stages. These are explored below.

Shelf water volume

Shelf water volume: A measure of the volume of water bounded inshore of the shelf-slope front. In this analysis, shelf water is defined as all water having salinity <34 psu. The position of the shelf-slope front varies inter-annually with the higher shelf water values indicating the front being pushed further towards the shelf break.

high shv: front pushed towards sbf low shv: front pushed inshore (more slope water on shelf)

Hypothesis: Higher recruitment success correlated with years of higher shelf water volume in spring/summer. These months months may be particularly important as that is when spawning is occurring and the position of the sbf may influence the position of larvae (away from spawning grounds).

Additional variables in this dataset are shelf water temperature and salinity which may also be indicative of habitat conditions.

# Shelf water volume
shlfvol <- read.csv(here::here('data/shelf_water_volume/ShelfWaterVolume_BSB_update.csv'))

# wrangling date info, converting doy to date and month 
yrs <- as.vector(nrow(shlfvol))
shlfvol$Year <- as.character(shlfvol$Year)
for (i in 1:nrow(shlfvol)){
  yrs[i] <- strsplit(shlfvol$Year, ".", fixed = TRUE)[[i]][1]
}
shlfvol$year <- yrs
shlfvol <- shlfvol %>% mutate(date_= as.Date(Year.Day-1, 
                                             origin=paste0(year, "-01-01")), 
                              month= strftime(date_, "%m"), 
                              day=strftime(date_,"%d"), 
                              year = as.integer(year),
                              month = as.numeric(month))  


# Create shw vol by month w/lag
df.lag = shlfvol %>%
  group_by(year,month) %>% 
  summarise(mean_t = mean(ShW.T),
            mean_s = mean(ShW.S),
            mean_v = mean(ShW.Vol)) %>% 
  mutate(mean_t_lag2 = lag(mean_t,2),
         mean_t_lag3 = lag(mean_t,3), 
         mean_t_lag6 = lag(mean_t,6),
         mean_s_lag2 = lag(mean_s,2),
         mean_s_lag3 = lag(mean_s,3), 
         mean_s_lag6 = lag(mean_s,6),
         mean_v_lag2 = lag(mean_v,2),
         mean_v_lag3 = lag(mean_v,3), 
         mean_v_lag6 = lag(mean_v,6))
# Join in-situ bottom temps w/assessment recruitment estimate
df.join = dplyr::full_join(recruit, df.lag, by = join_by(year)) %>%
  dplyr::select(year, month, recruit_est, mean_t, mean_s, mean_v, 
                mean_t_lag2, mean_s_lag2, mean_v_lag2,
                mean_t_lag3, mean_s_lag3, mean_v_lag3) %>%
  tidyr::drop_na()
# See what months have data
sort(unique(df.join$month))

## [1]  7  8  9 10 11

hist(df.join$month) # will group into spring/summer fall/winter categories

## Shelf water volume no lag
ggplot2::ggplot(df.join,
                aes(x=recruit_est, y=mean_v)) + 
  geom_point(color = 'black') +
  stat_smooth(method = "lm",
              formula = y ~ x,
              geom = "smooth") +
  xlab('Recruitment estimate') +
  ylab('Mean shelf water volume') +
  labs(title = 'Shelf water volume') +
  stat_cor(aes(label=..rr.label..)) +
  theme_bw()

## Shelf water temperature no lag
ggplot2::ggplot(df.join,
                aes(x=recruit_est, y=mean_t)) + 
  geom_point(color = 'black') +
  stat_smooth(method = "lm",
              formula = y ~ x,
              geom = "smooth") +
  xlab('Recruitment estimate') +
  ylab('Mean shelf water temperature') +
  labs(title = 'Shelf water temperature') +
  stat_cor(aes(label=..rr.label..)) +
  theme_bw()

## Shelf water salinity no lag
ggplot2::ggplot(df.join,
                aes(x=recruit_est, y=mean_s)) + 
  geom_point(color = 'black') +
  stat_smooth(method = "lm",
              formula = y ~ x,
              geom = "smooth") +
  xlab('Recruitment estimate') +
  ylab('Mean shelf water salinity') +
  labs(title = 'Shelf water salinity') +
  stat_cor(aes(label=..rr.label..)) +
  theme_bw()

With lags

2 Months

## Shelf water volume 2 month lag
ggplot2::ggplot(df.join,
                aes(x=recruit_est, y=mean_v_lag2)) + 
  geom_point(color = 'black') +
  stat_smooth(method = "lm",
              formula = y ~ x,
              geom = "smooth") +
  xlab('Recruitment estimate') +
  ylab('Mean shelf water volume') +
  labs(title = 'Shelf water volume - lag 2 months') +
  stat_cor(aes(label=..rr.label..)) +
  theme_bw()

## Shelf water temperature  2 month lag
ggplot2::ggplot(df.join,
                aes(x=recruit_est, y=mean_t_lag2)) + 
  geom_point(color = 'black') +
  stat_smooth(method = "lm",
              formula = y ~ x,
              geom = "smooth") +
  xlab('Recruitment estimate') +
  ylab('Mean shelf water temperature') +
  labs(title = 'Shelf water temperature - lag 2 months') +
  stat_cor(aes(label=..rr.label..)) +
  theme_bw()

## Shelf water salinity 2 month lag
ggplot2::ggplot(df.join,
                aes(x=recruit_est, y=mean_s_lag2)) + 
  geom_point(color = 'black') +
  stat_smooth(method = "lm",
              formula = y ~ x,
              geom = "smooth") +
  xlab('Recruitment estimate') +
  ylab('Mean shelf water salinity') +
  labs(title = 'Shelf water salinity - lag 2 months') +
  stat_cor(aes(label=..rr.label..)) +
  theme_bw()

3 Months

## Shelf water volume 3 month lag
ggplot2::ggplot(df.join,
                aes(x=recruit_est, y=mean_v_lag3)) + 
  geom_point(color = 'black') +
  stat_smooth(method = "lm",
              formula = y ~ x,
              geom = "smooth") +
  xlab('Recruitment estimate') +
  ylab('Mean shelf water volume') +
  labs(title = 'Shelf water volume - lag 3 months') +
  stat_cor(aes(label=..rr.label..)) +
  theme_bw()

## Shelf water temperature 3 month lag
ggplot2::ggplot(df.join,
                aes(x=recruit_est, y=mean_t_lag3)) + 
  geom_point(color = 'black') +
  stat_smooth(method = "lm",
              formula = y ~ x,
              geom = "smooth") +
  xlab('Recruitment estimate') +
  ylab('Mean shelf water temperature') +
  labs(title = 'Shelf water temperature - lag 3 months') +
  stat_cor(aes(label=..rr.label..)) +
  theme_bw()

## Shelf water salinity 3 month lag
ggplot2::ggplot(df.join,
                aes(x=recruit_est, y=mean_s_lag3)) + 
  geom_point(color = 'black') +
  stat_smooth(method = "lm",
              formula = y ~ x,
              geom = "smooth") +
  xlab('Recruitment estimate') +
  ylab('Mean shelf water salinity') +
  labs(title = 'Shelf water salinity - lag 3 months') +
  stat_cor(aes(label=..rr.label..)) +
  theme_bw()

Annual

# Create shw vol by year w/lag
df.lag = shlfvol %>%
  group_by(year) %>% 
  summarise(mean_t = mean(ShW.T),
            mean_s = mean(ShW.S),
            mean_v = mean(ShW.Vol)) %>% 
  mutate(mean_t_lag2 = lag(mean_t,2),
         mean_t_lag3 = lag(mean_t,3), 
         mean_t_lag6 = lag(mean_t,6),
         mean_s_lag2 = lag(mean_s,2),
         mean_s_lag3 = lag(mean_s,3), 
         mean_s_lag6 = lag(mean_s,6),
         mean_v_lag2 = lag(mean_v,2),
         mean_v_lag3 = lag(mean_v,3), 
         mean_v_lag6 = lag(mean_v,6))
# Join in-situ bottom temps w/assessment recruitment estimate
df.join = dplyr::full_join(recruit, df.lag, by = join_by(year)) %>%
  dplyr::select(year,recruit_est, mean_t, mean_s, mean_v, 
                mean_t_lag2, mean_s_lag2, mean_v_lag2,
                mean_t_lag3, mean_s_lag3, mean_v_lag3,
                mean_t_lag6, mean_s_lag6, mean_v_lag6) %>%
  tidyr::drop_na()

## Shelf water vol 

ggplot2::ggplot(df.join,
                aes(x=year, y=mean_v)) + 
  geom_point(color = 'black') +
  geom_line(color = 'black') +
  xlab('Year') +
  ylab('Mean shelf water volume') +
  labs(title = 'Shelf water volume') +
  theme_bw()

ggplot2::ggplot() + 
  geom_line(data = df.join, aes(x=year, y=mean_s), color = 'red') +
  geom_line(data = df.join,aes(x=year, y=mean_t*1), color = 'blue') +
  ylim(30.0,34.0) +
  scale_y_continuous(name = 'Sh.Water Salinity', 
                      sec.axis = sec_axis(~./1, name = 'Sh.Water Temperature')) + 
  xlab('Year') +
  labs(title = 'Shelf water salinity/temperature') +
  theme_bw()

ggplot2::ggplot() + 
  geom_line(data = df.join, aes(x=year, y=mean_s), color = 'red') +
  xlab('Year') +
  ylab('Mean shelf water salinity') +
  labs(title = 'Shelf water salinity') +
  theme_bw()

## Shelf water vol no lag

ggplot2::ggplot(df.join,
                aes(x=recruit_est, y=mean_v)) + 
  geom_point(color = 'black') +
  stat_smooth(method = "lm",
              formula = y ~ x,
              geom = "smooth") +
  xlab('Recruitment estimate') +
  ylab('Mean shelf water volume') +
  labs(title = 'Shelf water volume') +
  stat_cor(aes(label=..rr.label..)) +
  theme_bw()

## Shelf water temperature no lag
ggplot2::ggplot(df.join,
                aes(x=recruit_est, y=mean_t)) + 
  geom_point(color = 'black') +
  stat_smooth(method = "lm",
              formula = y ~ x,
              geom = "smooth") +
  xlab('Recruitment estimate') +
  ylab('Mean shelf water temperature') +
  labs(title = 'Shelf water temperature') +
  stat_cor(aes(label=..rr.label..)) +
  theme_bw()

## Shelf water salinity no lag
ggplot2::ggplot(df.join,
                aes(x=recruit_est, y=mean_s)) + 
  geom_point(color = 'black') +
  stat_smooth(method = "lm",
              formula = y ~ x,
              geom = "smooth") +
  xlab('Recruitment estimate') +
  ylab('Mean shelf water salinity') +
  labs(title = 'Shelf water salinity') +
  stat_cor(aes(label=..rr.label..)) +
  theme_bw()

With lags

## Shelf water vol 6 yr lag

ggplot2::ggplot(df.join,
                aes(x=recruit_est, y=mean_v_lag6)) + 
  geom_point(color = 'black') +
  stat_smooth(method = "lm",
              formula = y ~ x,
              geom = "smooth") +
  xlab('Recruitment estimate') +
  ylab('Mean shelf water volume') +
  labs(title = 'Shelf water volume - lag 6 years') +
  stat_cor(aes(label=..rr.label..)) +
  theme_bw()

## Shelf water temperature 6 yr lag
ggplot2::ggplot(df.join,
                aes(x=recruit_est, y=mean_t_lag6)) + 
  geom_point(color = 'black') +
  stat_smooth(method = "lm",
              formula = y ~ x,
              geom = "smooth") +
  xlab('Recruitment estimate') +
  ylab('Mean shelf water temperature') +
  labs(title = 'Shelf water temperature - lag 6 years') +
  stat_cor(aes(label=..rr.label..)) +
  theme_bw()

## Shelf water salinity 6 yr lag
ggplot2::ggplot(df.join,
                aes(x=recruit_est, y=mean_s_lag6)) + 
  geom_point(color = 'black') +
  stat_smooth(method = "lm",
              formula = y ~ x,
              geom = "smooth") +
  xlab('Recruitment estimate') +
  ylab('Mean shelf water salinity') +
  labs(title = 'Shelf water salinity - lag 6 years') +
  stat_cor(aes(label=..rr.label..)) +
  theme_bw()

Food availablity

Larval tilefish eat zooplankton, likely calanus. Calanus finmarchicus are a copepod (crustacean) with a one-year life cycle and are an important food source for many commercially important species. Calanus spp. are lipid rich, herbivorous species that eat phytoplankton, diatoms in particular (Hobbs et al. 2020).

Diatoms are often represented as microplankton (>20 µm), but many species are of the nanoplankton size class (2-20 µm), and a smaller few may even overlap with picoplanton size class (<2 µm).

Calanus

Calanus is not as common in MAB, need to figure out dominant zooplankton in MAB.

# Calanus
calanus <- ecodata::calanus_stage %>% filter(Time %in% c(1998:2021))%>% 
    rename_all(., .funs = tolower) %>% 
   mutate(year = time)


ggplot() +
  geom_line(data = calanus %>% filter(epu == 'GB', 
                                  var == 'adt-Spring'), 
       aes(x = year , y = value, col = epu), lwd = 1) + 
  geom_line(data = calanus %>% filter(epu == 'MAB', 
                                  var == 'adt-Spring'), 
       aes(x = year , y = value, col = epu), lwd = 1) +
  labs(color = c('EPU')) +
  theme_minimal()

Georges Bank

# GB
c5.gb <- calanus %>% filter(epu == 'GB', var == 'c5-Spring')
adult.gb <- calanus %>% filter(epu == 'GB', var == 'adt-Spring' )

df.c5 <- dplyr::full_join(recruit, c5.gb, by = join_by(year)) %>%
  dplyr::select(year, recruit_est, value) %>%
  tidyr::drop_na()
df.adt <- dplyr::full_join(recruit, adult.gb, by = join_by(year)) %>%
  dplyr::select(year, recruit_est, value) %>%
  tidyr::drop_na()

# Regression
ggscatter(df.c5, x = 'recruit_est', y = 'value', 
          add = "reg.line", conf.int = TRUE, 
          cor.coef = TRUE, cor.method = "pearson",
          xlab = "Recruitment estimate", 
          ylab = "Calanus c5 spring (No. per 100m^-3)",
           title = 'c5')  

ggscatter(df.adt, x = 'recruit_est', y = 'value', 
          add = "reg.line", conf.int = TRUE, 
          cor.coef = TRUE, cor.method = "pearson",
          xlab = "Recruitment estimate", 
          ylab = "Calanus adult spring (No. per 100m^-3)",
           title = 'Adult')  

# GLM 
eqn <- as.formula(paste('recruit_est ~', paste(colnames(df.c5)[1], 
                                               collapse = " + ")))

mod0 <- glm(recruit_est ~ 1, 
            data = df.c5, 
            family = "poisson")

mod1 <- glm(eqn, 
            data = df.c5, 
            family = "poisson")

summary(mod0)

## 
## Call:
## glm(formula = recruit_est ~ 1, family = "poisson", data = df.c5)
## 
## Deviance Residuals: 
##     Min       1Q   Median       3Q      Max  
## -1007.0   -634.1   -150.4    471.9   1112.8  
## 
## Coefficients:
##               Estimate Std. Error z value            Pr(>|z|)    
## (Intercept) 14.1979354  0.0001802   78775 <0.0000000000000002 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## (Dispersion parameter for poisson family taken to be 1)
## 
##     Null deviance: 8933474  on 20  degrees of freedom
## Residual deviance: 8933474  on 20  degrees of freedom
## AIC: 8933809
## 
## Number of Fisher Scoring iterations: 4

summary(mod1)

## 
## Call:
## glm(formula = eqn, family = "poisson", data = df.c5)
## 
## Deviance Residuals: 
##     Min       1Q   Median       3Q      Max  
## -1080.8   -584.0   -147.9    658.1    944.2  
## 
## Coefficients:
##                Estimate  Std. Error z value            Pr(>|z|)    
## (Intercept) 43.50429247  0.05804867   749.4 <0.0000000000000002 ***
## year        -0.01459584  0.00002891  -504.8 <0.0000000000000002 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## (Dispersion parameter for poisson family taken to be 1)
## 
##     Null deviance: 8933474  on 20  degrees of freedom
## Residual deviance: 8677368  on 19  degrees of freedom
## AIC: 8677705
## 
## Number of Fisher Scoring iterations: 4

AIC(mod0, mod1) %>% dplyr::arrange(AIC)

##      df     AIC
## mod1  2 8677705
## mod0  1 8933809

null_prediction <- exp(predict(mod0))
mod_prediction <- exp(predict(mod1))


plot(df.c5$year, df.c5$recruit_est, type = 'l')
lines(df.c5$year, null_prediction, col = "gray")
lines(df.c5$year, mod_prediction, col = "red")

Mid-atlantic

# Mid-Atlantic Bight
c5.mab <- calanus %>% filter(epu == 'MAB', var == 'c5-Spring')
adult.mab <- calanus %>% filter(epu == 'MAB', var == 'adt-Spring' )

df.c5 <- dplyr::full_join(recruit, c5.mab, by = join_by(year)) %>%
  dplyr::select(year, recruit_est, value) %>%
  tidyr::drop_na()
df.adt <- dplyr::full_join(recruit, adult.mab, by = join_by(year)) %>%
  dplyr::select(year, recruit_est, value) %>%
  tidyr::drop_na()

# Regression
ggscatter(df.c5, x = 'recruit_est', y = 'value', 
          add = "reg.line", conf.int = TRUE, 
          cor.coef = TRUE, cor.method = "pearson",
          xlab = "Recruitment estimate", 
          ylab = "Calanus c5 spring (No. per 100m^-3)",
           title = 'c5')  

ggscatter(df.adt, x = 'recruit_est', y = 'value', 
          add = "reg.line", conf.int = TRUE, 
          cor.coef = TRUE, cor.method = "pearson",
          xlab = "Recruitment estimate", 
          ylab = "Calanus adult spring (No. per 100m^-3)",
           title = 'Adult')  

Microplankton

# Microplankton

Chlorophyll-A

# CHL-A

References:

Hobbs, L., Banas, N. S., Cottier, F. R., Berge, J., & Daase, M. (2020). Eat or sleep: availability of winter prey explains mid-winter and spring activity in an Arctic Calanus population. Frontiers in Marine Science, 7, 541564.